High-energy magnon dispersion and multi-magnon continuum in the two-dimensional 

Heisenberg antiferromagnet 
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We use quantum Monte Carlo simulations and numerical analytic continuation to study high-energy 
spin excitations in the two-dimensional S = 1/2 Heisenberg antiferromagnet at low temperature. 
We present results for both the transverse (x) and longitudinal (z) dynamic spin structure factor 
Sa;,z(q, w) at q = (tt,0) and (ir/2,ir/2). Linear spin-wave theory predicts no dispersion on the line 
connecting these momenta. Our calculations show that in fact the magnon energy at (tt, 0) is 10% 
lower than at (n/2,n /2). We also discuss the transverse and longitudinal multi-magnon continua 
and their relevance to neutron scattering experiments. 

PACS: 75.40.Gb, 75.40.Mg, 75.10.Jm, 75.30.Ds 
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It is now well established that the spin-wave theory 
of the two-dimensional Heisenberg model jl] correctly 
describes the low-energy dynamics of layered antiferro- 
magnets such as La 2 Cu04 |§,[§. Several methods have 
been used to calculate the quantum renormalization fac- 
tor Z of the spin wave velocity, c = Zcq, with the result 
Z w 1.18 for S = 1/2 Early neutron scattering ex- 

periments were consistent with this renormalization over 
the entire Brilloin zone J9|. However, recently more accu- 
rate measurements have shown clear deviations at high 
energies in La2Cu04 jl(| and other materials (Til]. In 
particular, spin-wave theory predicts no dispersion on 
the magnetic zone boundary, i.e., along the line of mo- 
menta q = (jr — x, x), whereas the experimental data 
show a significant variation in the excitation energy. This 
could be caused by interactions other than the nearest- 
neighbor super-exchange J normally used to model the 
materials [[Io],[l2|. The deviations could also be at least 
partially due to a failure of low-order spin-wave theory 
to account for the high-energy dynamics of the Heisen- 
berg model. A series expansion calculation |ljj indeed 
indicates a momentum dependent renormalization fac- 
tor, with the magnon energy at q = (tt, 0) approximately 
7% lower than at (tt/2,tt/2). This, however, is opposite 
to the trend found experimentally for La2Cu04 pQ] . 

Another important issue, that has not yet been con- 
fronted experimentally, is the multi-magnon continuum 
expected to be present in the dynamic structure factor 
5(q, uj) above the single-magnon mode. The neutron 
scattering cross-section has been analyzed assuming only 
a delta-function representing the spin waves with disper- 
sion cj q P,p[-pT|. The resolution of present experiments 
is not sufficient for detecting the presence of additional 
spectral weight, much less to determine its distribution. 
Theoretical calculations have so far also not been success- 
ful in accurately determining the full dynamic structure 
factor. In order to provide guidance for improved fitting 
procedures it is imperative to obtain quantitative esti- 
mates of the multi-magnon spectral features, in particu- 



lar considering that experiments may soon become suf- 
ficiently accurate for detecting the continuum. Within 
spin-wave theory, correctly accounting for the interac- 
tions that give rise to the continuum is extremely com- 
plicated and has led to contradictory results P, |l^ , |l5[ . 
Numerical calculations of 5(q, to) are challenging because 
large lattices have to be used to converge to the thermo- 
dynamic limit, and the extraction of real-time dynamics 
from quantum Monte Carlo (QMC) data is difficult. Such 
calculations have therefore so far only given limited in- 
sights pq , p7[ . The series work Q suggests a significant 
continuum but gives no information on its shape. 

Here we address the issues of high-energy magnon dis- 
persion and multi-magnon continuum using QMC calcu- 
lations in a way that explicitly separates the transverse 
and longitudinal components of the dynamic structure 
factor. This allows us to more accurately determine both 
the magnon energy and the continuum part of the spec- 
trum. We consider the Heisenberg model on a square 
lattice, defined in standard notation by the Hamiltonian 



H = J^S,; 



(1) 



The coupling J « 1500 K in typical planar cuprates. 
On an infinite lattice the spin-rotational symmetry of 
this model is spontaneously broken at T — and long- 
range order develops along an axis that we take as the 
z-direction. The two-point correlations then have dif- 
ferent longitudinal (z) and transverse (x, y) components. 
In a finite system the symmetry is not broken and all 
correlation functions are equal to the same rotational av- 
erage. We want to access separately the transverse dy- 
namic structure factor, which contains the single-magnon 
mode as well as a continuum, and the longitudinal com- 
ponent comprising only a multi-magnon continuum. To 
this end, we explicitly break the symmetry by applying 
a staggered magnetic field; 



H 



(2) 



1 



We adjust the field strength h so that the induced stag- 
gered magnetization equals the known value in the ther- 
1(5?) | = 0.307 jiJ]|l9Uf-§. In the 
(at T = 0) and the 



modynamic limit; m 
limit of large system sizes h 



correlation functions become equal to their values in a 
system with spontaneously broken symmetry. 

In real layered cuprates there is a small inter-layer cou- 
pling and some degree of anisotropy, resulting in a finitc- 
T transition to an ordered state, in La2Cu04 at Tn « 300 
K. We here consider a low temperature, (3 — J/T = 32, 
corresponding to approximately 50 K, where the stag- 
gered magnetization in real materials is very close to the 
saturated T = value. Our staggered field h can then 
also be viewed as a mean-field treatment of the inter-layer 
coupling [^o| . 

Working with lattices with N = L x L spins and pe- 
riodic boundary conditions, we have used the stochastic 
series expansion method [^lj to calculate the imaginary- 
time dependent spin-spin correlation function 



G a (q,r) = (S a (~q,T)S a (q,0)), 



where a = x, z and 



S a (q,r) = * ]Te-^e-^fe^. (4) 



The correlation function is related to the dynamic struc- 
ture factor according to 



Ga(q, r) 



du>S a (q, uj)e 



(5) 



which in principle can be inverted to yield the real- 
frequency dynamics from the correlation function com- 
puted in the simulations. With simulation data affected 
by statistical fluctuations there are well known difficul- 
ties in carrying out this analytic continuation in practice, 
and one can only expect limited frequency resolution us- 
ing, e.g., the maximum entropy method ^2| . Here we will 
instead assume reasonable functional forms for the trans- 
verse (x) and longitudinal (z) dynamic structure factors, 
with parameters that are adjusted to satisfy the equality 
(||). With the QMC method used, derivatives of G a (q, r) 
can also be directly calculated (2^] and impose additional 
constraints on S a (q, u>); from Eq. (0) 
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We use the first (n = 1) and (n = 2) second derivatives, 
which at r = correspond to the first two frequency 
moments of the spectrum. We also use the sum rule 

, . 2 f°° duj „ . , 
X«(q) = - / — S a {q,u>), (7) 

TT J-oo U 

where Xa(q) is the static susceptibility 



x«(q) 



/ drG a {q, r). 
Jo 



(8) 



The spectrum also obeys the bosonic relation 
S a (q,-w)=e-^S a (q,u/). ' 

For the model forms of the transverse and longitudi- 
nal structure factors we take simple functions that reflect 
the expected gross spectral features. The transverse com- 
ponent should include a delta- function at an energy w qi 
representing the magnon, and a continuum that does not 
extend below w q (temperature broadening of the high- 
energy magnons should be insignificant at (3 = 32). The 
longitudinal component is entirely in the continuum and 
can also not extend below w q . Both continua must decay 
to zero rapidly as to — > oo, to ensure convergence of all 
frequency moments. We use 

S x (q,cj) = Ai(q)*(w-w fl )+j4 a (q)/»(q,£«>), (9) 
S z (q,u) = B(q)f z (q,Lj), (10) 



(3) where the continua f x ,z(q, are given by 



f x (q,Lu) = r x e-^-») ^ , (0 for uj < w q ), (11) 
f z {q,L0) = r z (uj ~ cugfe" 1 ^-^ , (0 for u < u q ), (12) 

where r x and r z are factors normalizing the frequency 
integrals to unity. For clarity, we have suppressed the 
dependence of the parameters v 1 a,p, a, and b on q. The 
static structure factors (the integrated spectral weights) 
are given by 

S x (q) - (S"(-q)S*(q)) - ^(q) + A 2 (q), (13) 
S z (q) = (S*(-q)S*(q)) = B(q). (14) 

In scattering experiments with unpolarized neutrons, a 
rotational average of the transverse and longitudinal dy- 
namic structure factors is measured. The rotationally av- 
eraged static structure factor 5(q) = 2S x (q) + S z (q) — 
2Ax + 2A 2 + B. 

The values we find for the "self-consistent" staggered 
field (demanding m = 0.3070 to an accuracy of better 
than 0.00005) are h/J = 0.10652 for L = 4, 0.022575 for 
L = 8, 0.005450 for L = 16, and 0.001615 for L = 32. 
With the QMC algorithm [^l] implemented in the z- 
basis, and with the field also in the z-direction, the lon- 
gitudinal correlations can be easily calculated. In order 
to compute the transverse correlations it is more efficient 
to apply the field in the cc-direction, so that the measure- 
ments can be carried out with diagonal operators [ pi[ . 
Hence, we perform two independent simulations for each 
lattice size. 

Fig. [l] shows our results for the ratio of the longitudinal 
and transverse spectral weights versus the momentum, 
along a standard path in the Brilloin zone. At q = (tt, tt) 
the ratio diverges (in the thermodynamic limit), reflect- 
ing the divergence of S z (tt, tt) due to the long-range order. 
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FIG. 1. Ratio of the longitudinal and transverse static 
structure factors (integrated spectral weights) along a path 
in the Brilloin zone. Open circles are for L — 16 and solid 
ones for L = 32. 



For q — > (0, 0) and q — ► (71", 71") it approaches zero and 
hence the magnon completely exhausts the total spectral 
weight in these limits (other treatments have shown that 
the magnon exhausts the low-energy transverse spectral 
weight ^5j). The ratio is the largest at q = (7r/2,7r/2) 
where it is above 0.7 — it is between 0.6 and 0.7 over 
much of the Brilloin zone. Hence, in unpolarized neutron 
scattering experiments, 30 — 35% of the cross-section at 
the zone boundary is due to the longitudinal continuum. 

Next we study the full dynamic structure factor, fo- 
cusing on the zone-boundary momenta q = (n, 0) and 
(7r/2,7r/2). The analytic continuation using fits to the 
functions (0) and ( |lo|) was carried out using QMC 
imaginary-time data with a spacing At = 0.5 up to r = 4 
(the data become too noisy for higher r). For L = 4, 
we have compared with exact diagonalization results and 
find that the magnon energy is very accurately repro- 
duced (better than 0.5%) but the weight in the magnon 
is underestimated by about 5% because the form of the 
continuum, Eq. (^), is not appropriate for a very small 
system where only a small number of significant delta- 
functions represent the continuum. Our results for larger 
systems should have smaller systematic errors. 

Fig. ^| shows the size dependence of the magnon weight 
and energy (the scatter of the data gives an indication of 
the statistical errors). The L = 4 data is from exact diag- 
onalization; in this case the two momenta are degenerate 
due to equivalence of this lattice and a 2 4 hyper-cube. 
Rough extrapolations in 1 / L give the infinite-size energy 
luJJ « 2.15 for q = (tt,0) and 2.39 for (tt/2, tt/2). The 
relative weight of the magnon in S x (q,uj) is w 60% at 
(tt, 0) and 85% at (tt/2, tt/2). Within linear spin-wave 
theory, including the quantum-renormalization Z = 1.18, 
the energy ui q /J = 2.36 for both momenta. Hence, the 
results show that the interactions neglected in spin- wave 
theory have strong effects at q = (tt,0), lowereing the 
energy by almost 10% and transferring « 40% of the 
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FIG. 2. Properties of the single-magnon part of the 
transverse dynamic structure factor vs inverse system size 
(L = 4,8,12,16,20,24, and 32). (a) single-magnon fraction 
of the weight, (b) magnon energy. 



magnon weight into the continuum. At (tt/2, 7r/2) the 
excitation energy is instead slightly increased by the in- 
teractions and a much smaller fraction of the weight is in 
the continuum. 

Fig. H shows the magnon as well as the transverse and 
longitudinal continua for L = 32. The rotationally aver- 
aged structure factor measured in neutron scattering ex- 
periments with unpolarized neutrons is also shown. The 
exponent p in the longitudinal continuum, Eq. (|l2|), is 
small and positive, p < 0.05, giving a spectral weight 
concentrated close to the magnon. The longitudinal con- 
tinuum is significantly narrower at (tt/2, tt/2) than at 
(tt, 0). For both momenta, the transverse continuum is 
broad and centered further above the magnon. All these 
spectral features are consistently present for system sizes 
L > 8. As seen in the insets of Fig. |[ for L = 32 only 
« 45% of the rotationally averaged spectral weight is 
in the magnon mode at (tt,0), increasing to w 65% at 
(tt/2, tt/2). The magnon weight decreases only slightly 
for even larger systems, as can be inferred from Fig. 

In previous work using QMC and analytic continua- 
tion [ fu||l~7| , the magnon energy was extracted as the first 
moment of the rotationally averaged relaxation function, 
F(q,uj) ~ S(q,Lj)/u}. Using this procedure, we find the 
moment 2.53 at (tt,0) and 2.58 at (tt/2, tt/2), i.e., both 
significantly above the actual magnon energies (the first 



3 



moment of S(c[, lo) is higher yet; approximately 2.66 for 
both momenta). This is clearly due to the large spectral 
weight in the continuum. There is potentially a simi- 
lar problem in the analysis of neutron scattering data. 
Since the present-day energy resolution is such that a 
non-negligible part of the continuum that we have found 
here would be included in an experimental fit to a single 
resolution-broadened peak, the extracted magnon energy 
may be too high. This effect would be particularly im- 
portant at and close to q = (it, 0), where our results show 
that less than 50% of the total weight is in the magnon 
peak and the continuum is relatively broad. 

Experimental results do appear to show the presence of 
some weight consistent with a continuum above the reso- 
lutioned broadened magnon peak but the statistics 
is not sufficient for extracting its size and shape. Our 
calculations should be useful for analyzing future exper- 
iments with higher resolution. Such experiments will be 
very important for determining the significance of other 
interactions in the antifcrromagnetic cuprates beyond the 
nearest-neighbor exchange J. The significant effects of 
magnon interactions that we have found at q = (tt, 0) 
may also have important consequences for the broaden- 
ing of the Raman spectrum Jl5| and the momentum space 
anisotropy of the photoemission spectrum |^6| . 
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FIG. 3. The transverse (solid curves) and longitudinal 
(dashed curves) dynamic structure factor at q = (n, 0) (upper 
panel) and (tt/2, n/2) (lower panel). S x (q,uj) and Sz(q, cj) are 
normalized to 1 and S z (q)/S x (q) [see Fig. pL respectively. 
The delta-function piece of S m (q, w) is represented by a verti- 
cal line with height (indicated by a star) equal to its relative 
weight A\ (q) . The insets show the rotationally averaged spec- 
tra, S(q,u>) = 2S a (q,o;) + S z (q,iu). 
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